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Abstract 

A method for including a priori information in the 2-D D-bar algorithm is presented. Two methods 
of assigning conductivity values to the prior are presented, each corresponding to a different scenario 
on applications. The method is tested on several numerical examples with and without noise and is 
demonstrated to be highly effective in improving the spatial resolution of the D-bar method. 


1 Introduction 

Electrical impedance tomography (EIT) is a low-cost, portable, and noninvasive imaging modality that is free 
of ionizing radiation with many potential applications for pulmonary imaging. In EIT, an image is formed 
of the conductivity distribution cr inside a body using only surface voltage and current measurements. 
Mathematically, this is a nonlinear inverse problem which is well known to be extremely ill-posed. A 
significant challenge in EIT imaging is the computation of static images with high-quality spatial resolution. 
Due to the ill-posedness, finer details in the image are often lost in the presence of noisy measurements. 
Including prior information in the reconstruction algorithm has been shown to be one way to improve spatial 
resolution This prior knowledge corresponds to a clinical situation in which 

we have a CT scan (or other similar data) for a human subject from which we may extract information 
regarding spatial locations of organ boundaries or conductivity estimates. When diagnosing and treating 
certain lung conditions, it is often necessary to obtain repeated thoracic CT scans, each of which imparts a 
dose of harmful radiation. EIT scans, on the other hand, have no ill effects. It is therefore highly desirable 
to use a priori information obtained from a CT or other scan to provide an improved EIT image, and then 
perform repeated harmless and comparatively inexpensive EIT scans in place of follow-up CT scans. 

Reconstruction algorithms that involve the minimization of a cost functional, such as a Gauss-Newton 
algorithm, include the a priori information in the penalty term, penalizing reconstructions that deviate too 
greatly from the prior in a given norm. This technique does not generalize to noniterative algorithms, and 
until now there has been no direct reconstruction method to utilize a priori information. The algorithm 
presented here therefore represents the first direct reconstruction method for EIT to incorporate a priori 
data. 

In this paper, we first provide an outline of the D-bar method without a priori data in © The a priori 
scheme for D-bar methods is described in SjJJ in which spatial information regarding locations of inclusion 
boundaries and approximated conductivity values are encoded into the equations for D-bar. Results using 
simulated data with and without noise on a circular domain with adjacent current patterns and 32 electrodes 
are presented in (j4j 

2 Theoretical background 

In EIT, current is applied on electrodes on the surface of a domain, and the resulting voltages are measured 
on the electrodes. Let H C K 2 be a bounded, simply connected Lipschitz domain. Then the electric potential 
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u(x , y ) within 12 is modeled by 


V • y)\7u(x,y)) = 0, (x,y)&Q, 


Cl) 


where a(x,y) is the conductivity distribution within 12. The boundary data for the inverse problem is given 
by the Dirichlet-to-Neumann (DN) map which takes the boundary voltages to the current densities on 
the boundary: 


A<t : u(x,y )| an ->• a(x,y 


an 


( 2 ) 


where v is the outward normal to the surface. The inverse conductivity problem is to reconstruct er(a;, y) for 
(x, y) £ fi given knowledge of the DN map A^. We will denote Ai to be the DN map corresponding to the 
case of constant conductivity er = 1. 

In this work, we will modify the 2-D D-bar reconstruction method based on the constructive global 
uniqueness proof in [21] to include a priori information about the conductivity. In m it is shown that twice- 
differentiable conductivities er can be uniquely determined from A a . It is well-known that for er £ C 2 (f2), 
equation (JT]) can be transformed to the Schrodinger equation through the change of variables q = 
u = \fo. Under the assumption that er is constant in a neighborhood of cA2, one can smoothly extend 
q = 0 to M 2 and consider the Schrodinger equation in the entire plane. Points (x, y) £ M 2 will be identified 
with points z = x + iy £ C. The direct reconstruction method relies on exponentially growing solutions 
to this Schrodinger equation, also known as complex geometrical optics (CGO) solutions, which arise when 
a complex frequency parameter k = k\ + ik 2 £ C is introduced in the equation and one seeks solutions 
ip(z,k) satisfying the asymptotic condition e~ lkz tp(z,k) -1 6 tU 1,p (]R 2 ),p > 2. It was proved in |2T] that 
the Schrodinger equation 

(—A + q(z))ip(z,k) = 0, zeM 2 (3) 


has a unique solution with the desired asymptotic property when q is of the form 

The CGO solution if)(z, k ) and its relative fi{z, k) := e~ lkz ip{z, k ) are key to the direct reconstruction 
algorithm. The algorithm has been implemented numerically and tested on simulated and experimental data 
EiiEuisniHiiianra. a nonlinear regularization method for this algorithm was provided with proof 
in m- Since details of the algorithm appear in numerous places in the literature, such as the papers above 
and |19j . we will give only a brief summary here with a focus on computation and the regularized method. 
As in m we will assume without loss of generality in this section that er = 1 in a neighborhood of dfl. 

The steps of the algorithm are 

Step 1. Compute the CGO solution ip on cA2 from the DN data by solving the boundary integral equation 



\ d n = e. lkz \dn- f G k (z - C)(A ff - Ai)^-, k)ds, 
J an 

where G k is the Faddeev Green’s function m for the Laplacian, defined by 

G k (z) := e ikz [ -A G k (z) = S(z). 

Step 2. Compute the scattering transform from 

t(fc) = f e lkz (A a — Ai )ip(z, k ) ds. 

Jan 

Step 3. Solve the D-bar (5) equation for p,(z, k ) 

d k F(z,k) = 

47 rk 

which can be written in integral form as 


1 

(2^F 



t(fc') 

k'(k — k') 


e- i(zfc ' +2/? V(z,fc/) dk’. 


(4) 

(5) 

( 6 ) 

(7) 


F(z, k) = 1 + 


( 8 ) 
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Step 4. Compute the conductivity <j(z) for each z £ 0 or in a region of interest from 

a(z) = p 2 (z, 0), z £ Q. (9) 

In the regularized method, the functions i/j(z,k) in Step 1, t(k ) in Step 2, and p(z 7 k) in Step 3 are 
computed for complex frequencies |fc| < R , where R serves as a regularization parameter J7j dependent 
upon the noise level. In an ideal setting with no noise present, the reconstruction converges to the actual 
conductivity as R —> oo pointwise E3- Our a priori method will take advantage of this fact, balancing the 
ability to choose R large for ideal priors with fidelity to the data. 


3 Outline of the a priori method 


To motivate the a priori scheme, we first note that the scattering transform can be written in terms of 
a scattering transform computed from a prior known conductivity distribution er pr in the form t (k) = 
t pr + perturbation. To this end, given <7 pr , let A pr denote the DN map corresponding to er pr , let tp pl denote 
the CGO solution satisfying 


tV( 2 , A:)Ian = e lfc2 |pn - [ G k {z - C)( A P r - Ai)^ pr (-, k)ds(Q, 
Jan 


let t pr denote the scattering transform satisfying 


tpr(fe) = / e lkz (A pi - Ai)^ pr (z, k)ds(z), 


ia n 


and let /i pr denote the CGO solution satisfying 

These equations are valid for cr pr £ L°°(fl) by (2- Subtracting (111 from § , we see 


( 10 ) 


( 11 ) 


( 12 ) 


t(fc)-t pr (fc) = 


e (A pip A i'll) A pr i/> pr F Ai^ pr )ds 


i an 


[ e^ 2 (A a (l/j - V’pr) - Ai(V> - Ippr) + (A ct - A pr )V> pr )ds. 
Jan 


Thus, 


t(fc) = tj 


r (k) + ( t 
Jan 


e lkz [(A a - Ai ){ip - ^ pr ) + (A a - A pr )^ pr }ds. 


(13) 


Formula (13) suggests the following scheme. Given er pr , compute A pr from a numerical forward solver, such 
as FEM, compute tp pr and t pr from (10) and 0, respectively, compute t (k) from ( [l3| , and use this t (k) 
in the D-bar method. However, this natural approach has several drawbacks when applied to noisy data. 
First of all, since the measured data A^ has noise, it is necessary to compute t (k) on a truncated domain 
\k\ < R. This means finer details encoded in large |fc| values of the prior will be lost. Second, the numerical 
computation of A pr itself introduces error that is not necessarily a good match to the noise in A a . Thus, the 
term 

f e z ' ks (A a -A pi )^ pr ds 
Jan 

is not an accurate perturbation, and errors in A a — A pr will be amplified by the exponentially growing 
functions e lkz and tp pl . 

An alternative approach motivated by (13) is to define an approximation to the scattering transform 
piecewise by 

ft(fc), \k\<Rx 

Ir 1 ,r 2 (^) := S tpr (k), R\ < |fc| < 1?2 ■ (14) 

l o, |fc| >122 
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In this approximation, the perturbation term in (13) is neglected for |fc| > R\, and the entire scattering 
transform is truncated for some R 2 > R\. Neglecting this term for |fc| > R\ is motivated by the fact that the 
size of R\ is limited since t(fc) will inevitably blow up for larger values of |fc| in the presence of noisy data. 
Since the computation of t pr is noise-free and in general much more numerically robust than the computation 
of t, we may select R 2 to be significantly larger than R\. The larger the value of f? 2 , the stronger the influence 
of the a priori information. The next question is then how to compute t BliH2 . For |fc| < f?i, t(fc) can be 
computed using which is the same as computing t EIE , using the notation in m- To avoid the problems 
that arise from computing A pr , the scattering transform t pr can be computed directly from the definition 
of the scattering transform, provided cr pr £ C 2 . Then, defining q pl = —^=-, the scattering transform t pr is 
defined to be the nonlinear Fourier transform of q pl {2T| 


t pr (k) := / e llcz q pl (z)ip pr (z,k)dz, 


(15) 


where ip pl is the solution of the Schrodinger equation with q = q pr . Once the scattering transform has 
been computed, the CGO solution p can be solved from (|8|. We define p R2 as the solution to 


p R2 (z, k) = ! + 


lJ R 1 ,R 2 


(fc') 


(27t) 2 J\k\<R 2 h?{k-k?) 


e -i(«fc+ 2 fcO A ( 2!)fc / ) dk \ 


(16) 


However, there is one more thing to note. The Green’s function for the D-bar operator dk is Ar, and so the 
solution ^ to ([7]) is obtained from 


p(z, k) = lim 


1 


R -¥00 I irR 2 




p(z, k)dk - 




l\k\ <R 




(17) 


where the first term tends to 1 as R —> 00 . Thus, 


i^R 2 (z,k) 


1 


1 


b R 1 ,R 2 


(A'') 


7Tl? 2 


/ p{z,k)dk , . - 

l\k\<R 2 (27r)^ J\k\<R 2 k f (k — k') 


i(zk'+zk') 


F R2 (z,k') dk 1 . 


(18) 


In practical computations, since p R2 (z, k) is unknown, the first integral is replaced by 1, that is, its limit as 
R 2 —^ 00 . Approximating p R2 in this term by /x pr , we have derived an equation for the approximation 


Hr 2 ( z , k ) = 


1 


1 


7Tl? 2 


/ Ppr(z, k)dk , . - 

'\k\<R 2 J\k\<R 2 k {k — k') 


t fll .fl 2 (fc ) e -i {z k'+zk' )llR ^ z k ,^ dk , ( 19 ) 


If the prior coincides with the correct conductivity distribution, this method converges and introduces no 
artifacts as R \, i? 2 —^► 00 by the convergence proof for the regularized D-bar method m- 

The strength of the prior, or its influence on the reconstruction, depends on i? 2 - If R 2 = -Ri, then 
t Kl , H , (k) = t® IE (fc), and the only influence of the prior on the reconstruction is in the term 


Mint (z) := 


1 


■kR 2 


' \k\<R 2 


M P r (z, k)dk. 


( 20 ) 


We can exert control over the amount of influence this term has by introducing a weighting parameter 
a £ [0,1] and writing 


Mr 2 ,«(z, k) = a + (1 — a)pi n t(z) + —^ / 

{zk) j |i 


e- i( ~ zk ' +Ep )p R2ta ( Z , k') dk', 


\k\<R 2 k'{k - k') 


( 21 ) 


which is equivalent to (16) if a = 1 and (19) if a = 0. 
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3.1 Computational considerations 


We now describe the numerical details for the computation of the conductivity distribution <t R 2 , 0 correspond¬ 
ing to the CGO solution /i R2ta to (21), including numerical approximations for the necessary intermediate 
operators and functions. 

As described in m , we compute a finite-dimensional matrix approximation L a to the DN map A a by 
first computing the discrete ND map R ff and then forming its inverse: L a = R^ 1 . Denote the number 
of linearly independent current patterns by N and the number of electrodes by L. To compute R<j, first 
orthonormalize the matrix of bipolar current patterns to obtain the set { JJ 71 }, m = 1,..., N and l = 1,..., L, 
and then apply the appropriate change-of-basis formula to the voltages, yielding {V™}. Then the ND matrix 
can be approximated by 


Re, (to, n) 


L 

E 

i=i 


^jrvr = Ev t j, 


where A is the area of an electrode and A 6 is the angular distance between electrodes. 

The discrete matrix approximation Li of the DN map Ai corresponding to homogeneous conductivity is 
obtained by first numerically solving the forward conductivity problem, using FEM for example, to create 
simulated voltage data for the case where a = 1 in fh The method described above for the computation of 
L „ may then be used to compute Li from this simulated data. 

The CGO solution is found by numerically solving the boundary integral equation 0 at the center 
zi of each electrode, as described in m- As in [I], we express the Faddeev Green’s function Gk as 


G k (z) 


—Re(EI(ikz)), 


where EI(z) is the exponential integral function, which in Matlab can be computed easily using the built-in 
function: EI(^) « 2 * EXPINT(z). To form a matrix approximation Tk for G k {zi — £);'), we must be careful 
of the logarthmic singularity that occurs when V = l. We therefore discretize the surface of each electrode 
into S points Zi s , s = 1,..., S, and compute 


r a n = J ^ Re ( EXPINT (* fc (~; - 0'))), i' ± I 

k ’ 1 2¥(El) Ef=i Rfi(EXPINT(iA:(2, - 0.))) I'= l ' 

Denote by bk = (&i(fc),..., bpf(k)) T and Ck the vectors of coefficients for the functions ip{z,k)\dn and 
e* fez |on, respectively, expanded in the basis of orthonormalized current patterns, so that if>{z,k)\an ~ 
Jbk,e lkz \an « Jck- We may then approximate the convolution of Gk with (A^ — A i)i/j for each z = Zi 
as a finite-dimensional vector: 


[ G k (z - Q(A a - AOOO, k)ds( C) « r k J(L CT - Li)b k , 
Jon 


/ an 

and the discrete version of Q is 

Jbk — Jck DkJ(L CT Li)bk- 

Multipying through by the transpose of the orthonormal matrix J yields the linear system 

[I + J T r k J(L (T -L 1 )]bk = c k , 


( 22 ) 


where I is the TV x N identity matrix. In our implementation, this system was solved in Matlab using the 
MLDIVIDE function. 

The computation of the CGO solution tp pl corresponding to the prior is handled quite differently. From 
m , we know /i pr satisfies the Lippmann-Schwinger type equation 


Mpr — 1 9k * (OprMpr)- 


(23) 


where gk(z) = e~ lkz Gk(z) is a fundamental solution of the operator —A — likd. The numerical solution of 
( |23| is based on ideas presented in [25] • and a complete description of the computational steps can be found 
in |19j. In short, we may write (23) as the linear system 


[I + g k * ? P r(-)]/V = 1, 


( 24 ) 
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which is solved for /i pr using a matrix-free method such as GMRES [2T or BICGSTAB [5B], which was used in 
our implementation, separating real and imaginary parts as required for the linear solver. The action of the 


linear operator on the left-hand side of (241 may be approximated efficiently using FFT and IFFT operations. 
Once we obtain the computation of t pr from its definition (15) is accomplished using simple numerical 


quadrature over the mesh of 2 -values. The integral /q nt is likewise found by applying numerical quadrature 
to the integrand p pr . 


computational methods described in [T9] . We write (21) as the linear system 


To obtain p R2a (z), we must solve the equation (21) for each z, which involves modification of the 


[I - AT{-)]n R2ta (z) = a + (1 - a)nint(z) 


(25) 


where the actions of the operators T and A are defined by 

t Hl ,R 2 (fc) j( z k'+zk 


Tf(k ) = 


Ink 


')/(fc), Ag{k) = - 


9(k) 


n J\k\<R 2 k ~ k' 


dk'. 


Note that (25) is not complex-linear due to the presence of the conjugate operator, so it is necessary to 


solve real and imaginary parts separately. This system is solved by again using a matrix-free solver such as 
BICGSTAB, where the action of A can be approximated by FFT and IFFT operations. 

From /J, R2 , a we obtain the resulting conductivity distribution a R2tC , = li R2ta (z, 0), which is a result of both 
the EIT data and the a priori data encoded into and /q nt . 


3.2 Constructing the prior: Blind Estimate Method 

We now discuss the first of two possible methods for assigning conductivity values to the a priori conductivity 
distribution <r pr . Each of these two methods describes the construction of a discontinuous a priori distribution 
<7 pr . To satisfy the requirement that <r pr £ C 2 (fl), we later mollify tj pr to obtain cr pr . 

In both methods, knowledge is assumed of the spatial locations of boundaries for various domain inclusions 
(such as heart, lungs, etc. for thoracic imaging) in the plane of the electrodes. In a clinical setting, this could 
be obtained by extracting the organ boundaries from a CT scan to obtain polygonal approximations to 
the actual organ boundaries. In our simulations we created polygonal boundaries representing heart, lungs, 
aorta, and spine within a circular domain, as shown in Figure [l] 



Figure 1: Simulated organ boundaries representing heart, lungs, aorta, and spine within a circular domain, 
used as the simulated a priori information in our experiments. Organ boundaries are approximated by 
polygonal regions; the dots represent polygon vertices. 
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In the blind estimate method for assigning conductivity values to the prior, we simply make educated 
guesses for the conductivity values within each approximate organ boundary. These values can, for example, 
be estimated from literature sources wherein conductivity values for human tissue have been reported. Let 
P C denote the polygonal region inside a particular approximate organ boundary and let {z r } be the finite 
set of points in the 2 -mesh used to construct cr pr . An approximate constant conductivity value o> is selected 
for P, and we assign <7 pr (z„) = o> for all z n G {z r } fi P. We repeat this process for all organ boundaries 
used in the prior, obtaining the conductivity distribution cr pr . Refinements to this process could be made by 
specifying regions of differing conductivities within individual organs if known inhomogeneities exist. 

The blind estimate method is much computationally simpler and faster than the alternative extraction 
method which will be described next. However, if any pathologies have developed between the time of the 
initial CT scan and the time of the EIT scan, these pathologies will not be reflected in the prior, and their 
expression in the final reconstruction is therefore entirely dependent on the EIT data. The full a priori 
scheme with the blind estimate method for constructing the prior is outlined in Algorithm [lj 


ALGORITHM 1 A priori scheme, Blind Estimate Method 
1: Obtain a priori information and EIT data: 

2: - Form polygonal approximations to organ boundaries. 

3: - Collect EIT data and compute L CT . 

4: - Use FEM to simulate homogeneous data and compute Li. 

5: Form computational grids for the k and z planes. 

6: Make blind estimates for conductivity values to form o- pr (z). 

7: Compute cr R2iCX : 

8: - Select Ri,R ,2 with Ri < R 2 . 

9: - Compute q pr = Ay/a/y/a. 

10: for \k\ < Ri do 

11: - Solve \22\ for to get ip\ qq ~ Jb^ 

12: - Compute t (k) from e-. 

13: end for 

14: for Ri < \ k\ < R 2 do 

15: - Solve ( |24| for /z pr = e~' lkz ip pT . 

16: - Compute t pr (fc) from d- 

17: end for 

18: for z E do 

19: - Solve ( |25| > for pR 2ta (z, •). 

20: - Compute cr R2:0i (z) from {§. 

21: end for 


3.3 Constructing the prior: Extraction Method 

In the extraction method, we first compute a reconstruction a from the EIT data alone using the D-bar 
method described in © Note that the first-order approximation t Gxp to the scattering transform as given 
in |14| could be used here as well to obtain an initial reconstruction. We then extract conductivity values 
from this reconstruction to obtain estimated values for the prior. This method has advantages over the 
blind estimate method in that pathologies not present in the CT scan data but that are apparent in the 
reconstruction a may be included in the prior. The full a priori scheme with the extraction method for 
constructing the prior is outlined in Algorithm [2j 

In our experiments with simulated data, we developed and used the following techniques for extracting 
approximate conductivity values from the reconstruction a to be assigned to the lungs, heart, aorta, spine, 
and background in the prior. In what follows, we denote by {z' s } the finite set of points in the 2 -mesh used 
to construct er, and the polygonal regions inside specific organ boundaries by P heart , P apine , etc. 


1. Lungs. We examine each lung in the a reconstruction and compare the appearance of the lungs to the 
a priori approximate lung boundaries. If, based on the reconstruction <r, the lungs appear to be free 
of pathology (i.e. there are no suspicious inhomogeneities within the regions Pi [ung and P rlung ), then 
the following method may be used. Find the set P) lung := {z' s } fl P l lung , and compute the average value 


M 


' 1 lung 


: = M E < Z 'm 


G Pn 


(26) 
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where M denotes the number of points z' m in P, i ung , and then assign er pr (~n) = cr, lung for all z n £ 
{z r }nP 1 lung . This process is then repeated for the right lung. On the other hand, if the reconstruction 
a reveals possible lung pathologies in the form of inhomogeneities within a lung region, then the 
method can be revised in the following way. Assuming (without loss of generality) that one or more 
inhomogeneities appear in the left lung, divide the region P, lung into a finite number of connected 
subsets Sj C Pi i ung , where each subset represents an area of fairly homogeneous conductivity in the 
reconstruction a. Then for each Sj, compute the average conductivity over the points z' m £ {z' s } fl Sj 
and assign this value to a pr (z n ) for z n £ {z r } D Sj. 

2. Heart and aorta. To compute cr pr values for the heart region, one could potentially employ a similar 
method to that described for the lung regions. However, the position and shape of the reconstructed 
heart is more sensitive to noise level and truncation radius than the lung, and therefore using the 
anatomical position in the prior may include extraneous pixels. See Figure [3] as an example of how the 
size and shape of the reconstructed heart can vary. The aorta, on the other hand, is typically invisible 
in the reconstruction a, so such a method could not be used to assign cf pr values within the aorta. The 
following method is therefore given as an alternative to the method used for the lungs. First, define 
the quantities 

<j max := max {<r( z' m )}, cr min := min { o(z ' m )}, 

z , rn ^:{z , s ]r\^i z' rn ^:{z , s }r[^i 

and compute the value r = cr min + c(cr max — a min ) where c £ (0.5,1) is selected empirically. A good 
choice for c should optimally result in the set H := {z' m £ {z' s } fl fl : cr(z' m ) > t} being selected so as 
to be a connected subset of fl and roughly the same size as the region P heart , and such a c may vary 
depending on noise levels and choice of truncation radius in the computation of a. We find the set H 
and compute 

1 M / 

^he a rt ■ 1 1 ^ ^ ® (An) ) An f 

m—1 

where M denotes the number of points in H. Finally, assign tf pr ( 2 n ) = cr hGart for all z n £ {z r } fl P hoart , 
and, since the aorta likely has conductivity values very similar to those of the heart, further assign 
these same values to cr pr in the region P aorta . 

3. Spine. Due to its small size, the reconstruction of the spine typically has very poor spatial resolution 
and its appearance and associated conductivity values can vary widely in the reconstruction a in the 
presence of noise. Since we can usually assume that the spine is one of the most resistive objects in a 
thoracic EIT scan, we simply assign cr pr (Ai) = cr min for all z n £ {z r } n P spina . 

4. Background. We define the background to be the set Pb g fl — U jPj where each P 7 corresponds to 
an organ boundary included in the prior, and assign values to P bg according to the following method. 
Compute the quantities n = <r min + Ci(cr max - cr min ) and r 2 = cr mio +c 2 (cr max - <r min ) where Ci,c 2 £ (0,1), 
Ci < c 2 . Find the set B := {z' m £ {z'| : cr(z' m ) £ [t!,t 2 ]}, and compute 

1 M 

°bg ■= Z 'm £ B > 

m—1 

where M denotes the number of points in B. Assign a pr (z n ) = cr bg for all z r £ { z n } fl P bg - As with the 
value c used for the heart, the values Ci and c 2 must be selected empirically, and may once again vary 
depending on noise levels and choice of truncation radius in the computation of cr. Since the lungs, 
which have low conductivity compared to the background, tend to dominate the reconstruction, it is 
usually advantageous to choose Ci,c 2 to be skewed to the upper end of the scale (0,1). Well-chosen ci 
and c 2 should result in the set B excluding most of the region corresponding to the lungs and spine, 
as well as the high conductivity region inside the heart. 
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ALGORITHM 2 A priori scheme with Extraction Method 
1: Obtain a priori spatial information and EIT data: 

2: - Form polygonal approximations to organ boundaries. 

3: - Collect EIT data and compute L^. 

4: - Use FEM to simulate homogeneous data and compute Li. 

5: Compute conductivity o-(z) using standard D-bar methods: 

6: - Form computational grids for the k and 2 : planes for both a and cr pr 

7: - Select a truncation radius R±. 

8 : for |fc| < Ri do 

9: - Solve \22\ for b k to get ijj\dn « Jb k 

10: - Compute t (k) from §■ 

11: end for 

12: for 2 : G O do 

13: - Solve the (Ri-truncated) equation § for p(z, •). 

14: - Compute < 7 ( 2 :) from (j9j. 

15: end for 

16: Extract conductivity values from cr(z) to form cr pr . 

17: Compute (Tfl 2jQ (z): 

18: - Compute q pr = Ay'o’/y'o’. 

19: - Select R 2 > R±. 

20: for Ri < |/c| < R 2 do 

21: - Solve ( |23| > for ipp r (-,k). 

22: - Compute t pr (fc) from < |15| ). 

23: end for 

24: -FormtH 1 ,R 2 - 

25: - Select a. 

26: for 2 : G £2 do 

27: - Solve ( |25| ) for pR 2 , a {z, •). 

28: - Compute (Jr 2 , ci .{z) from § 

29: end for 


3.4 Iterative approaches 


The a priori schemes described in the previous pages may be used alone to obtain a reconstruction 
but there is potential for further refinement of spatial resolution through the use of iterative approaches. 
The motivation is to take advantage of the enhanced spatial resolution in a R2ta to construct a new prior that 
is more accurate than the original in terms of conductivity values and possible pathologies. We may include 
in this updated prior any new information that appears in the reconstruction <t R 2jQ . This may be especially 
advantageous if the blind estimate method was used to construct the original prior, but the patient has since 
developed some pathology that is visible in the EIT data. Another situation where iteration may provide 
enhanced results is if we desire to use the extraction method, but the reconstruction a has very poor spatial 
resolution. In fj4]we provide an example of the first of these scenarios, using simulated data. 

The computational steps in these iterative approaches are the following: (1) obtain the reconstruction 
<7 H2iQ , using either the blind estimate or extraction method to assign conductivity values to the prior, (2) use 
the extraction method described in 13.3 to extract conductivity values from o R a (rather than from a), (3) 
use these extracted conductivity values to form an updated prior a' vr , (4) repeat the a priori scheme using 
the original EIT data with the updated prior o' to obtain an updated reconstruction o' R a . This entire 


process could potentially be repeated again if desired to obtain a second iterate o" R 


4 Results 

We now present the results of our test problem using simulated data. In this test problem, we simulate 
a situation in which a priori information is available from a previous CT scan, but the patient has since 
developed a pleural effusion in one lung. We tested the previously described a priori scheme using both the 
blind estimate and extraction methods for assigning conductivity values, and an iteration step as described 
above was performed on the results from the reconstructions using the blind estimate method. 

We assumed that the organ boundaries without pleural effusion were given by the polygonal approxi¬ 
mations shown in Figure [l] and we assumed a circular domain of radius 143.2 mm. We therefore created 
a phantom with these same organ boundaries, domain shape, and dimensions. Conductivity values were 
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assigned to the phantom heart, lungs, aorta, and spine, and the FEM method with the complete electrode 
model including contact impedance described in [2D] was used to generate voltage data. To simulate a pleural 
effusion, conductivity was increased in the phantom in the bottom of the left lung. We will use the conven¬ 
tion that the left lung appears on the left-hand side of the image. The phantom with assigned conductivity 
values is shown in Figure [2] Random zero-mean Gaussian noise was added to the simulated voltages at 
0%, 0.1%, and 0.2% of the maximum voltage values; the D-bar reconstructions a using the method of ^2] 
for each of the three noise cases are shown in Figure [3] All reconstructions in this section, including the 
reconstructions a and the results of the a priori schemes, were computed using a z-mesh with 101 x 101 
elements, R\ = 3.8, and we tested values R 2 £ {3.8,5.0,7.5,10}, and a £ {0,0.5,0.75,0.9}. Note that the 
a priori organ boundaries are the correct boundaries, and the conductivity values and distributions will be 
modified. 



Figure 2: Phantom representing pleural effusion. Conductivity values are in S/m. 
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(a) Noise level = 0%. 


(b) Noise level = 0.1%. 


(c) Noise level = 0.2%. 


Figure 3: Plots of the reconstructions cr computed using the regularized D-bar method of jj2] (see also lines 
[5][l5]of Algorithmic with truncation radius R\ = 3.8, at noise levels 0%, 0.1%, and 0.2%, with superimposed 
actual organ boundaries. Each noise case is plotted on its own scale; these scalings will be used for all plots 
concerning each noise level within this paper. 



4.1 Blind estimate method applied to test problem 

We assigned “blind estimate” a priori conductivity values representing a phantom with two homogeneous 
lungs with conductivity 0.200 S/m, in contrast to the actual values displayed in Figure [C The values for the 
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background, heart, aorta, and spine differed slightly from the actual values. These “blind estimates” were 
used for all three noise cases, and are given in Tabic [lj Reconstructions using the blind estimate method 
can be seen in Figures |6j [9j [l2j 

Given the obvious lung pathology apparent in the reconstructions, we then performed an iteration step 
wherein the left lung was divided into two regions, which we shall refer to as the “lung top” and “lung 
bottom,” separated by a horizontal line segment. We computed conductivity values for the iterate a' pl . 
separately in each of these two regions, using the methods described in §3.3| to extract conductivity values 
from the ctr 2 , q reconstruction with R 2 = 5.0 and a = 0.75. In the computation of the values for the heart, 
aorta, and background, we selected values c = 0.85, ci = 0.25, and C 2 = 0.95; for simplicity, we used these 
same values in all noise cases. The resulting conductivity values used in <7p r are also given in Table [lj 

The location of the dividing line between lung top and lung bottom was chosen by visually inspecting 
the cr ii2iQ reconstructions and selecting a horizontal line at which to form the division. For simplicity, we 
used the same approximate dividing line in all three noise cases. This approximate dividing line is compared 
to the actual lung division used to create the phantom in Figure [4] The reconstructions resulting from the 
iteration step can be seen in Figures [TJ [lOj and [13} 



Figure 4: Locations of dividing line between the “lung top” and “lung bottom.” The dividing line used in the 
phantom is indicated by a solid line. The approximate dividing line (dashed line) was used in the extraction 
method and the iteration step for the blind estimate method, and was obtained by visually inspecting the 
u R2 ,oi reconstructions from the blind estimate method. 


It is evident from the reconstructions in Figures [6[ |9[ [T2| that with or without noise the blind estimate 
method without iteration detects the pleural effusion provided a very small value of a is not combined with 
a small value of i? 2 - The case a = 0 and R 2 = 3.8 corresponds to using only fi- mt as in equation (191 and 
the scattering transform from the regularized D-bar method without a prior. Increasing a and decreasing 
R 2 weakens the influence of the prior. In the case of a strong prior, the prior dominates the reconstruction 
in the blind estimate method, resulting in good spatial resolution of organ boundaries, but poor detection 
of the effusion. Adding the iteration step described in §3.4| results in excellent detection of the effusion in 
every case, and the aorta and spine can be clearly seen with excellent spatial accuracy even with a prior 
of medium weight, such as a = 0.5, R 2 = 7.5. The ringing effect seen in the reconstructed spine and lungs 
for small a or large R 2 is likely due to the influence of the term /ij nt on the reconstruction. This term 
provides excellent spatial resolution of the prior, detecting edges extremely well, but introducing ringing 
since it becomes increasingly uniform, tending to 1 as i ?2 ~^ 00 . This effect can be seen in Figure [5j 


4.2 Extraction method applied to test problem 

In assigning approximate conductivity values, for each noise level we first reconstructed a with R\ = 3.8 
using the D-bar method with no a priori information. The reconstruction of a is plotted along with the a 
priori organ boundaries in Figure [3] 

For the extraction of conductivity values, from the a reconstructions with the superimposed a priori 
organ boundaries, it was clear that the left lung has increased conductivity toward the bottom, so we again 
divided the lung into top and bottom to construct the prior. For simplicity, we used the same dividing line 
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(a) R 2 = 3.8 (b) R 2 = 5.0 (c) R 2 = 7.5 (d) R 2 = 10 

Figure 5: Plots of the real part of /q nt used in the simulations with various truncation radii R 2 . Since 
Mint —> 1 as i ?2 -» oo, the scale must be adjusted for each value of R 2 for best viewing results. 


as was used in §4.1[ and we selected the same values for c, Ci, and C 2 for all noise cases. Using the methods 
outlined in |3.3[ we extracted conductivity values from the reconstruction a to create er pr ; these assigned 
values are given in Table [2j along with the conductivity values used in the phantom for comparison. We then 
proceeded with the rest of the a priori scheme outlined in fj3j using i?i = 3.8, and testing various values for 
R 2 and a. The resulting reconstructions are given in Figures [Hj [TT] and 

In this method, with or without noise, the effusion is clearly visible for all weights of the prior, with 
improvement in the organ shapes as the weight of the prior increases. Excellent reconstructions are found 
even in the presence of noise. An iteration step is not included for this method since the first step produces 
very high quality reconstructions. 


14 


5 Conclusion 

A method for including a priori information in the 2-D D-bar algorithm was presented with two methods 
suggested for assigning conductivity values to the prior. The a priori information is included in the scattering 
transform and in the integral equation for the CGO solution /i(z, k) and is weighted with two parameters R 2 
and a in the scattering transform computation and the integral equation for p, respectively. The method 
is demonstrated to be highly effective on numerically simulated data with noise levels typically used in EIT 
data simulations. The method shows promise for clinical use in lung imaging when a priori information 
about organ boundaries can be obtained from a recent CT or MRI scan, for example. Future work is needed 
to evaluate its clinical effectiveness. 

Table 1: Conductivity values in S/m for the phantom as well as the “blind estimate” er pr values assigned, 
along with values assigned to <r pr in the subsequent iteration step, for each of the 3 noise cases. 



Back¬ 

ground 

Heart 

L Lung 
top 

L Lung 
bottom 

R Lung 

Aorta 

Spine 

Used in phantom 

0.424 

0.750 

0.240 

0.600 

0.240 

0.750 

0.150 

“Blind estimates” used in cr pr 

0.500 

0.800 

0.200 

0.200 

0.200 

0.800 

0.100 

Used in <Tp r , 0% noise 

0.431 

0.798 

0.261 

0.364 

0.233 

0.798 

0.187 

Used in cr pr , 0.1% noise 

0.427 

0.767 

0.274 

0.333 

0.233 

0.767 

0.178 

Used in <jp r , 0.2% noise 

0.450 

0.858 

0.247 

0.428 

0.247 

0.858 

0.163 


Table 2: Conductivity values in S/m for the phantom as well as cx pr values computed using extraction method 
for each of the 3 noise cases. 



Back¬ 

ground 

Heart 

L Lung 
top 

L Lung 
bottom 

R Lung 

Aorta 

Spine 

Used in phantom 

0.424 

0.750 

0.240 

0.600 

0.240 

0.750 

0.150 

Extracted from cr, 0% noise 

0.401 

0.681 

0.283 

0.398 

0.251 

0.681 

0.186 

Extracted from cr, 0.1% noise 

0.393 

0.648 

0.292 

0.373 

0.252 

0.648 

0.178 

Extracted from cr, 0.2% noise 

0.423 

0.742 

0.272 

0.460 

0.260 

0.742 

0.177 
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Blind Estimate Method 
0% Noise 
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<0 
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No Prior, R t = 3.8 


Strong Prior 


Weak Prior 





q = 0.9, Ri =10 q = 0.75, Ri = 10 a = 0.5, R 2 = 10 a = 0, ft, = 10 



0.1 0.2 03 0.4 0.5 06 0.7 08 09 



a = 0.9. R 2 = 5.0 a = 0.75, R 2 = 5.0 a = 0.5, R t = 5.0 n = 0. R t = 5.0 



Figure 6: Reconstructions <J R2ta for the 0% noise case using the a priori scheme with the blind estimate 
method (before iteration step), with various values of a and R 2 . The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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Blind Estimate Method 
Plus Iteration Step 
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a = 0.9. Ri = 3.8 a = 0.75, R 2 = 3.8 a = 0.5. R 2 = 3.8 o = 0. R 2 = 3.8 



a = 0.9, R 2 =10 a = 0.75, R 2 = 10 a = 0.5, R 2 = 10 a = 0. R 2 = 10 
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Figure 7: Reconstructions <j' R2 a for the 0% noise case using the a priori scheme with the blind estimate 
method plus one iteration step, with various values of a and R 2 - The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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Extraction Method 
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Figure 8: Reconstructions ct H 2iQ for the 0% noise case using the a priori scheme with the extraction method, 
with various values of a and f? 2 • The reconstruction with no prior is at the top for comparison. The strength 
of the prior increases moving left to right and top to bottom. The scale bar at the bottom applies to all 
reconstructions at this noise level. 
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Blind Estimate Method 
0.1% Noise 
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Figure 9: Reconstructions (J R2 , a for the 0.1% noise case using the a priori scheme with the blind estimate 
method (before iteration step), with various values of a and R 2 - The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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Blind Estimate Method 
Plus Iteration Step 
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a = 0.9. Rj = 3.8 a = 0.75, R 2 = 3.8 r» = 0.5. R 2 = 3.8 o = 0. R 2 = 3.8 



a = 0.9, R 2 =10 q = 0.75, R 2 = 10 a = 0.5, R 2 = 10 q = 0, R 2 = 10 
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Figure 10: Reconstructions a' Rct a for the 0.1% noise case using the a priori scheme with the blind estimate 
method plus one iteration step, with various values of a and R 2 - The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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No Prior, R i = 3.8 


Extraction Method 
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Figure 11: Reconstructions cr H2i0 , for the 0.1% noise case using the a priori scheme with the extraction 
method, with various values of a and 1 ? 2 - The reconstruction with no prior is at the top for comparison. 
The strength of the prior increases moving left to right and top to bottom. The scale bar at the bottom 
applies to all reconstructions at this noise level. 
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Blind Estimate Method 
0.2% Noise 
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a = 0.9, Rj = 3.8 a = 0.75, R 2 = 3.8 a = 0.5. R 2 = 3.8 a = 0. R 2 = 3.8 





a = 0.9, R 2 =10 a = 0.75, R 2 = 10 a = 0.5, R 2 = 10 a = 0, R 2 = 10 



Figure 12: Reconstructions cr R2 0L for the 0.2% noise case using the a priori scheme with the blind estimate 
method (before iteration step), with various values of a and R 2 - The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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Blind Estimate Method 
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Figure 13: Reconstructions a' Rct a for the 0.2% noise case using the a priori scheme with the blind estimate 
method plus one iteration step, with various values of a and i? 2 - The reconstruction with no prior is at the 
top for comparison. The strength of the prior increases moving left to right and top to bottom. The scale 
bar at the bottom applies to all reconstructions at this noise level. 
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Extraction Method 
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Figure 14: Reconstructions <jR 2 , a for the 0.2% noise case using the a priori scheme with the extraction 
method, with various values of a and R 2 . The reconstruction with no prior is at the top for comparison. 
The strength of the prior increases moving left to right and top to bottom. The scale bar at the bottom 
applies to all reconstructions at this noise level. 
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